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Abstract 

In a recent paper (El Omari and Le Guer, IJHMT, 53, 2010) we have investigated mixing and 
heat transfer enhancement in a mixer composed of two circular rods maintained vertically 
in a cylindrical tank. The rods and tank can rotate around their revolution axes while their 
surfaces were maintained at a constant temperature. In the present study we investigate the 
differences in the thermal mixing process arising from the utilization of a constant heat flux as 
a boundary condition. The study concerns a highly viscous fluid with a high Prandtl number 
Pr = 10, 000 for which this chaotic mixer is suitable. Chaotic flows are obtained by imposing 
temporal modulations of the rotational velocities of the walls. By solving numerically the 
flow and energy equations, we studied the effects of different stirring protocols and flow 
configurations on the efficiency of mixing and heat transfer. For this purpose, we used 
different statistical indicators as tools to characterize the evolution of the fluid temperature 
and its homogenization. Fundamental differences have been reported between these two 
modes of heating or cooling: while the mixing with an imposed temperature results in 
a homogeneous temperature field, with a fixed heat flux we observe a constant difference 
between the maximal and minimal temperatures that establish in the fluid; the extent of 
this difference is governed by the efficiency of the mixing protocol. 
Keywords: Convection, Imposed wall heat flux, Chaotic mixing, Unstructured finite 



* Corresponding author, Tel: +33 559 407 147, Fax: +33 559 407 160 
Email addresses: kamal.elomari@univ-pau.fr (Kamal El Omari), yves.leguer@univ-pau.fr (Yves 
Le Guer) 



Preprint submitted to International Journal of Heat and Mass Transfer 



December 14, 2010 



volume method. 



Nomenclature 

C p heat capacity (J/kg.K) 

k thermal conductivity (W/m.K) 

L wall characteristic length (m) 

p pressure (Pa) 

q surface heat flux (W/m 2 ) 

R 3 tank radius (m) 

Ri,R,2 rod radii (m) 

t time (s) 

T temperature (i(T) 

C/ tangential velocity (m/s) 

U velocity field (m/s) 

Dimensionless numbers 

Nu Nusselt number 

Pe Peclet number 

Pr Prandtl number 

Re Reynolds number 

X rescaled dimensionless temperature 

Greek symbols 

e eccentricity (m) 

p fluid density (kg/m 3 ) 

a standard deviation 

t period of modulation (s) 

r viscous stress tensor (Pa) 

Q angular velocity (rad/s) 

Subscript 

c cell 

m mean 

/ face of a cell 

2 



1. Introduction 

Mixing processes are widely encountered in many practical engineering domains where 
enhancements of the heat, mass and momentum transfer are required. In the present study, 
mixing was achieved through the existence of a laminar chaotic flow, which ensures the 
efficient stretching and folding of the material lines (here for a two-dimensional (2D) flow). 
The need for chaotic mixing is particularly interesting when high viscosity and/or shear- 
sensitive fluids are concerned [1]. In this case, classical laminar 2D time- independent flows 
are unable to give a good mixing performance and as a consequence, heat transport from the 
wall will be ineffective [2]. Such a situation exists for molten polymers or polymer blends [3], 
in food processing [4] or are encountered for the formation of highly concentrated emulsions 
[5] . Also, the processing of highly viscous polymer molecule networks or biological fluids in 
classical turbulent flow mixers can cause the degradation of these fluids in the high-shear 
regions of the flow. Besides, undesired or ill-defined fluid structures may be obtained. 

Chaotic flows were successfully used for heat transfer [6], especially in twisted-pipe heat 
exchangers [7] . Even if the flow in these pipes is steady, the chaotic character is obtained by 
elbows succession, where at sufficiently high Dean numbers, secondary vortices appear and 
promote chaotic advection [8, 9]. This type of mixer is efficient at intermediate Reynolds 
numbers, roughly for 40 < Re < 2000, and is thus not suitable for highly viscous fluids such 
as the one addressed in the present study. 

The thermal chaotic mixing has some particularities compared to the mixing of concen- 
tration (dye for example). Even if concentration and temperature are both passive scalars, 
their mixing require different strategies. The origin of this difference is the location of the 
scalar source. In the case of mixing of concentration, the two constituents are present in the 
mixer at the start of the process. Whereas, for the thermal mixing, the fluid is initially at a 
uniform temperature (a cold temperature for example in the case of a heating process) and 
the heat source is located at the mixer boundaries (in the absence of volumetric heating as 
the Ohm's effect). Two different situations are then encountered depending on the type of 
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the boundary condition imposed at the mixer walls: (i) if a fixed temperature is imposed, 
the amount of energy extracted from these walls by the fluid (parietal heat flux rate) depends 
on the ability of the flow to generate fluid trajectories perpendicularly to them and to mix 
this hot fluid with the colder one far from the boundaries, (ii) if a fixed flux is imposed at 
the walls, the rate of heat flux is then known and the behavior of the flow will impact the 
repartition of this heat between the different zones of the fluid. However, in both situations, 
the flow in the mixer must insure a good transport from the walls towards the center. 

In the case of concentration mixing no source is present since the second constituent is 
initially injected in the mixer as a small blob in the case of dye mixing for example, or it 
occupies a larger portion of the mixer space in the case of the blend of two fluids. As a 
consequence, the mixing strategies classically developed for concentration mixing may not 
be successful for thermal mixing if it neglect the mixing from and near the walls. 

When mixing a blob of dye in an enclosed fluid, Gouillart et al. [10] showed that the 
presence of the no-slip wall of a non-rotating tank slows down the rate of mixing, while the 
rate is exponential when this tank is rotating [11]. This is due to the presence of parabolic 
points in the non-rotating wall which brings unmixed fluid to their zone of interest, i.e. 
the center of their mixer. When the tank is rotating, the central zone is isolated from the 
peripheral unmixed zone and the rate of the scalar decay is exponential as expected for a 
fully chaotic flow. When dealing with thermal mixing as it is the case in the present work, 
the situation is reversed since the scalar (temperature) needs to be driven from the wall. 
We showed in [12] that the presence of parabolic points on the boundaries is a necessary 
condition in order to mix fluid from the walls and bring it towards the mixer center. 

The mixer geometry used here was derived from the classical eccentric cylinder geometry 
(Journal bearing flow) considered by different authors [13, 14, 15, 16]. This geometry pro- 
duces a large recirculation bubble, which is difficult to avoided even with a careful choice of 
stirring protocols. In contrast, the two-rod mixer studied here is suitable for obtaining a full 
chaotic flow without KAM regions. It is similar to the geometry of the two-roll-mills studied 
in the literature by Price et al. [17] or Young et al. [18] and Chiu et al. [19]. This geometry 
was studied in details by Jana et al. [20] under the name of "Vortex mixing flow" . They used 
numerical and experimental tools to study the mixing of a scalar (dye in the experiments) 
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by rotating alternately the two rods in a fixed tank (the outer cylinder). In our case, the 
tank has the ability to rotate. We previously studied this mixer in the case of heating of a 
Newtonian fluid using a fixed temperature at the mixer walls [21, 12]. We extended these 
studies to the case of shear-thinning and shear-thickening fluids [22]. 

The objective of the present work is to analyze the combined effects of mixing and heat 
transfer for the particular wall boundary condition of constant heat flux (i.e. the Neumann 
condition). In the authors' previous studies [21, 12, 22], mixing and heat transfer were 
studied for the constant wall temperature boundary condition (i.e. the Dirichlet condition) 
in the same mixer. In such a case, the imposed wall temperature represents an asymptotic 
limit for the evolution of the mean fluid temperature, and the rapidity to reach this limit is 
controlled by the efficiency of the mixing. In contrast, for constant wall heat flux, there is 
no asymptotic limit for the evolution of the mean fluid temperature because its evolution is 
prescribed by the imposed parietal heat flux density. As a consequence, the stirring strategy 
and the choice of mixing indicators must be selected in agreement with the type of wall 
heating considered. The differences in mixing induced by the choice of boundary condition 
(constant wall heat flux versus constant wall temperature) will be highlighted throughout 
this work. 

2. Problem formulation 

2.1. Geometry of the two-rod rotating mixer 

A sketch of the mixer studied is presented in Fig. 1. It is composed of two circular rods 
of equal radii, which are maintained vertically inside a cylindrical tank (a bounded domain). 
The tank and the rods are heated and can rotate around their respective revolution axes. 
The rods and the tank rotate alternately or with a continuous modulation of velocities. 
They can also have different directions of rotation. We have shown for a Newtonian fluid 
[21, 12] that this two-rod mixer is suitable for obtaining full chaotic flow without KAM 
regions, which is particularly interesting for industrial applications. The geometry of the 
flow domain is characterized by the radii of the rod and the cylindrical tank, which are 
respectively: R\ — R2 — 10 mm and -R3 = 50 mm. The eccentricity of the rods is set to 
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the value of e = 25 mm. The choice of e = -R3/2, i.e. the rods located at the midpoint 
between the center of the tank and its boundary, was found to be close to the optimum in 
our previous study [12] . This value is also near to the optimum of eccentricity found by Jana 
et al. [20] for a radii ratio R3/R1 = 4 and for a fixed tank and alternately rotating rods. 




Figure 1: Sketch of the two-rod mixer. 

2.2. Fluid properties 

The Newtonian fluid considered in this study has the thermophysical properties listed in 
Table 1. 



Dynamic viscosity (/x) 


1.5 Pa.s 


Density (p) 


990 kg.m~ 3 


Thermal conductivity (A) 


0.15 W.m~ l K~ l 


Specific heat (C p ) 


1000 J.kg^K" 1 


Prandtl number (Pr) 


10 4 


Peclet number (Pe) 


16,584 



Table 1: Fluid properties. 



Considering these properties, the thermal diffusivity a is equal to 1.515 x 10 7 m 2 .s 1 . 
The maximum angular wall velocities are fixed to: fii = = 30 rpm for the rods and 
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r2 3 = 6 rpm for the outer tank. The tangential wall velocity is then the same, and equal to 
U = 31.41 mm.s -1 . A characteristic Reynolds number for this flow can be evaluated as: 

In our study, the Peclet number Pe is large (see Tab. 1) hence the need to accelerate the 
mixing of the temperature by advection is evident. 

2.3. Flow parameters: stirring protocols 

Chaotic mixing flows are produced by varying the angular velocity of the rods and tank 
in time using a sine-squared waveform. The stirring protocols studied are defined by two 
parameters: the respective direction of the rotation between the rods and tank, and the 
size of the period of the velocity modulation. Other parameters, such as the ratio between 
the maximum velocities of the rods and tank, and the amplitude of the angular velocity 
modulation, are not considered here. The three possible flow configurations giving rise to 
different flow topologies are specified in Table 2. The rods are always rotating in phase. The 
symbol (+) indicates a counter-clockwise direction and (— ), a clockwise direction. 



Flow configuration 


Rod 1 


Rod 2 


Tank 


A 


(+) 


(-) 


(+) 


B 


(+) 


(+) 


(+) 


C 


(") 


(-) 


(+) 



Tabic 2: The three possible flow configurations related to the sign of angular velocity. 

In order to facilitate comparisons, the same types of temporal modulation are considered 
as in our previous work [12]. These are, for chaotic flows: a continuous modulation denoted 
CM (a sine-squared variation of the wall velocities) and a non-continuous modulation (al- 
ternated) denoted ALT, i.e. the rods are stopped together for half a period while the outer 
tank is rotating; for the next half-period the contrary applies. Three periods of modulation 
t are considered here : 15 s, 30 s and 60 s. These stirring protocols are illustrated in Fig. 
2. The maximum angular velocity is the same for the two types of modulation. The rods 
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Figure 2: The temporal modulation of the angular velocities of the rods and tank for the continuous (CM) 
and non-continuous (ALT) stirring cases. Here the modulation period, r, is 30 s. 

and the tank follow sine-squared modulations, they are detailed in [12]. In our modeling, we 
take into account inertial effects as the flow considered for this tangential wall velocity does 
not satisfy the quasi-steady hypothesis. 

2.4- The constant wall heat flux boundary condition 

The Neumann boundary condition (constant wall heat flux) is imposed on the rods and 
tank walls. As mentioned above, in previous studies [21, 12], we considered the constant wall 
temperature boundary condition (i.e. Dirichlet condition). This new boundary condition 
(Neumann condition) imposes to reconsider the objective of the mixing and heat transfer 
optimization problem (i.e. the optimal choice of the stirring protocols). Indeed, while in the 
case of a constant wall temperature the final mean temperature is asymptotically bounded 
by the parietal temperature, in the case of a constant wall heat flux the mean temperature 
T m will instead increase linearly with time and be directly proportional to the parietal flux. 
Effectively, for a 2D mixer we have a slope of: 




(2) 



where V is the volume of the fluid and S the outer surface of the mixer. 

Thus, the statistical values to follow in order to characterize the efficiency of the process 
are not the same ones used previously. Instead of following the mean temperature, we follow 
the temporal evolutions of the difference between the maximum and minimum temperatures 
or the variation of one of them compared to the mean temperature (see section 5). There- 
fore, in the case of poor mixing, extremum temperatures (T max and T m j n ) should have high 
values in comparison with the case of a good mixing process. Higher temperatures will be 
encountered in the vicinity of the heating walls and may cause undesirable damages to the 
fluid. 



3. Computational modeling 

3.1. The governing equations 

The unsteady Navier-Stokes equations governing the flow of incompressible fluids and 
the continuity equation are considered in their integral form: 

^- [ pU dV + [ pUU -ndS = [ -Vj> dV + / r • ft dS (3) 
<n Jv Js Jv Js 



U-ndS = (4) 
where r, the viscous stress tensor for a Newtonian fluid: 



s 



t = p(yu + (yu) T ) (5) 

The integration is over a volume V surrounded by a surface S oriented by an outward 
unit normal vector n. The equation of the energy conservation is considered in terms of the 
temperature: 

^ I p C P T dV + [ pC p TU-ndS = ! kVT-ndS (6) 
ut Jv Js Js 
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where k is the thermal conductivity. The temperature is non-dimensionalized by: 

T* = T A ^ m with AT = kllJi ( 7 ) 
k 

AT is fixed to 1, thus, q re f = — with L = 2 (R 3 — Ri). 

LJ 

At the beginning of the mixing process, the fluid is at a uniform initial temperature T n j 
and the heat flux at the different walls is fixed to q — 10 q re f. 



3.2. The numerical method 



The conservation equations (3, 4 and 6) are solved by the means of an in-house code called 
Tamaris. This code has a three-dimensional unstructured finite volume framework that is 
applied to hybrid meshes. It is a direct extension to 3D of the 2D code used in our previous 
study [12]. The cell shapes can be of different forms (tetrahedral, hexahedral, prismatic 
or pyramidal). The variable values (U, p and T) are stored at cell centers in a collocated 
arrangement. This three-dimensional (3D) code can also deal with 2D computations (e.g. 
in (x, y) plan) without any change, by considering a single layer of computational cells (in 
z direction) and by neutralizing the top and bottom faces (with respect to z). In the scope 
of this work, all of the computational meshes were generated using the open-source software 
Gmsh [23] . More details about the numerical method and the discetization procedures used 
are reported in the appendix A. The code validation was undertaken in [22] and extended 
here in the appendix B with a supplementary experimental test case. 



4. Mixing indicators 



In order to quantify the efficiency of the heating process for all the fluid inside the mixer, 
we used the following instantaneous measure: the standard deviation a of the dimensionless 
temperature T*. This quantity is defined as (the summation is made over all the mesh cells 
c of area A c ): 

i 

2 

(8) 



a 



E (Mt: - T; n ?) 
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where in the non-dimensional mean temperature calculated as: 



The standard deviation is an indicator of the level of homogenization of the scalar tem- 
perature inside the 2D tank [12]. The evolution of the mean temperature gives the total 
energy supplied to the fluid during the mixing process, which is directly proportional to the 
value of the wall heat flux when we consider the Neumann boundary condition (see section 
2.4). 

Another mixing indicator is the temperature scalar dissipation indicator, which is defined 



as: 



^(x>ii*ne) 



X 9 = t-/ HVTTdS = ^— > AJ|VT*||c (10) 



Stot J Slot 



where Stot is the total surface of the fluid in the mixer. This indicator illustrates the 
mechanism of scalar gradient production and destruction in the fluid. 

Heat exchange at a wall (rod or tank) is characterized here by the mean Nusselt number 
(we recall that q — 10 q re f, see section 3.1): 

^=9^/ k(r L -T*) dS = J" I T* '-T^ 3 (U) 

Jwall J S wa u t{ '\ 1 wall 1 m) >~>wall J S waU 1 wall 1 rn 



5. Results and discussion 



The mesh of the two- rod mixer used in this work contains 10680 computational cells 
and uses regular quadrilateral cells near the walls to enhance the resolution of the boundary 
layers (the figure of this mesh was shown in [22] so it is not reproduced here). This mesh 
was adopted after a grid size-dependence study as detailed in [22] which was extended to 
the case of prescribed heat flux. 

In this section we describe the obtained temperature field for the three stirring protocols 
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studied and we show the impact of these protocols on the defined mixing indicators, tem- 
perature extrema and distributions (PDF). A particular attention is paid to the differences 
with the case of constant temperature boundary condition studied in [12]. 



5.1. Temperature patterns 

V 

0.00 0.163 0.325 0.488 0.650 




Figure 3: Dimcnsionless temperature fields at time t = 120 s for the flow configuration C and for: (a) 
non- modulated (NM), (b) continuously modulated (CM) and (c) alternated (ALT) stirring protocols (period 
of modulation r = 30 s). 

In Fig. 3, the dimensionless temperature fields are shown for the stirring condition C 
of Table 2 and for three cases of the modulation of the wall angular velocity. After four 
periods (t = 120 s), we find that the temperature fields of the non- modulated (NM) and 
continuously modulated (CM) cases are quite similar. They follow the shape of the roughly 
steady-state streamline patterns and present large zones of unmixed cold fluid. Due to 
the continuous rotation of the walls in these two cases, there are permanent streamlines 
parallel to the boundaries that inhibit heat transport by advection to the center of the 
mixer. The flow patterns in the two-rod mixer were exhaustively presented and discussed in 
our previous paper [12] dealing with imposed temperature as a boundary condition, so they 
are not reproduced here. Unlike the NM and CM cases, the alternated modulation (ALT) 
case allows hot fluid strips to penetrate towards the center of the mixer. Indeed, when one 
of the walls is at rest, stagnation points form at its boundary and give rise to such a hot 
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fluid streams. As said earlier, when the walls are non-rotating, parabolic points appear in 
them, and by alternated rotation they take place successively in the tank and rods surfaces. 
In Figure 3(c), we can indentify two hot streams originating from the tank boundary, carried 
by the rotation of the rods (at t = 120 s, the rods are rotating at their maximum velocity 
while the tank is at rest, see Fig. 2). Also, with an alternated rotation, quite different flow 
topologies with crossing streamlines succeed each other in the mixer, and this feature is the 
key of efficient chaotic mixing [1]. 

When comparing the temperature fields obtained for the fixed temperature boundary 
condition to those of Fig. 3, we find exactly the same patterns, but with different scales. 
For the former case, the maximum value is 1 and is always present in the field, while the 
minimum value evolves continuously towards 1. For the fixed flux boundary condition, both 
extremums are evolving. These observations and the impact of the mixing efficiency on them 
are detailed in section 5.4. 

5.2. Influence of the wall boundary condition on the temperature homogenization 



0.45 

0.4 
0.35 

0.3 
0.25 

0.2 
0.15 

0.1 
0.05 





...»•••**********"*"*" 

/ 


* 

• ■ 












J / 


ALT ----- 

CM 
NM 

T m 



500 



1000 
time [s] 



1500 




- 0.: 



0.6 



0.4 



- 0.2 



10 
2000 



1000 
time [s] 



2000 



(b) 



Figure 4: The temporal evolution of the temperature standard deviation for the non- modulated (NM), 
continuously modulated (CM) and alternated (ALT) stirring protocols for a modulation period of 30 s and 
the flow configuration C of Table 2. (a) The imposed heat flux, (b) The imposed temperature (note that in 
this case the a axis is in the logarithmic scale). 



The temporal evolutions of the standard deviation a of the temperature scalar and its 
mean value are shown in Fig. 4 for the three stirring protocols, NM, CM and ALT. The 
results of the imposed heat flux as a boundary condition are compared to those of the 
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imposed temperature for a modulation period of 30 s and for the same flow configuration C 
of Table 2. 

In the first instants (t < 100 s), the evolutions of a are the same because during this 
time, the mechanism of heat transfer is only governed by thermal diffusion near the wall 
boundaries. Afterwards, the evolutions of a differ according to the stirring condition and/or 
protocol chosen. 

In Fig. 4(a) (imposed heat flux), each curve of a displays a rapid exponential increase 
before reaching a plateau value. At the beginning of the mixing process, the temperature 
fluctuations in the mixer increase by much more rate than the increase in the mean value T^, 
resulting in an increase of a, and later, an equilibrium is reached (among stretching, folding 
and diffusion), indicating that the total variation T* — does not evolve any more over the 
mixer section (see also Fig. 7(a) hereafter). This equilibrium is reached more rapidly for the 
ALT protocol. This is a major difference when compared to the case of the constant wall 
temperature boundary condition (Fig. 4(b)), where we found that a decreases exponentially, 
indicating a tendency towards complete homogenization. As described earlier, while the 
evolution of in the first case is a straight line with a predefined slope (see section 2.4), 
the curves obtained in the second case present an asymptotic behavior towards the maximum 
value 1. The rate of this evolution depends on the efficiency of the mixing protocol. 

In Fig. 4(a) we observe that the lower a plateau value is obtained for the ALT stirring 
protocol, whereas the highest is obtained for the NM stirring protocol. Also 95% of the final 
value of u is obtained after 250s for the ALT stirring protocol against 1500s for the NM 
stirring protocol which represents a speed of homogenization of temperature fluctuations 
roughly 5 times faster. An asymptotic thermal regime is then more rapidly achieved for 
the ALT stirring protocol. This confirms what was observed in Fig. 3. For the constant 
wall temperature boundary condition (Fig. 4(b)) we have obtained the same result, i.e., the 
rate of the decrease of a is higher for the ALT stirring. We infer that in order to avoid 
the persistence of closed streamlines near the walls, it is necessary to alternately move these 
walls. Thus, the heat transport from the walls can be enhanced. Similarly to the NM stirring 
protocols, the CM stirring protocols also give rise to closed streamlines in the vicinity of the 
walls and these closed streamlines prevent the radial transport of the temperature scalar 
inside the mixer. For the NM and CM cases, close to the rotating boundaries, the radial 
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heat transfer is made only by conduction across the streamlines. 



5.3. Influence of flow configurations and of period length 



0.35 




time [s] 



Figure 5: Temporal evolution of temperature standard deviation for the continuously modulated (CM) and 
alternated (ALT) stirring protocols for a modulation period of 30 s and flow configurations A, B and C of 
Table 2. 

Mixing in the two-rod mixer can be conducted in different manners with respect to the 
direction of rotation of the three walls, as listed in Table 2. In Fig. 5, the evolutions of a 
are compared for the cases of the CM and ALT stirring protocols with a modulation period 
of 30 s and flow configurations A, B and C. Regardless of the flow configuration, the ALT 
stirring protocol always gives the best result (i.e., lower value of a). For this stirring protocol, 
the values of a are almost insensitive to the choice of the flow configuration. This behavior 
is also observed for the constant wall boundary condition (see Fig. 13 of reference [21]). 
The CM stirring protocol results, however, depend on the flow configuration, and the best 
efficiency is obtained with configuration A. 

Another parameter that impacts this mixing process is the length of the modulation 
period r of the angular velocity. Fig. 6 shows the influence of the duration of the modulation 
period on the a value for the ALT stirring protocol and flow configuration C. For this flow 
configuration, longer duration of the modulation period is related to lower a value as well 
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as better mixing and heat transfer efficiency. Small periodic oscillations are observed on the 
curves with the same period as the modulation period. Their amplitudes are highest for the 
longest modulation period (i.e. 60 s). An exhaustive study of the impact of the period size 
was reported in [12] (for Dirichlet boundary condition), which corroborates these findings 
and indicates that there is an optimum of r at about 60 s for configuration C and at about 
30 s for configuration B. 



0.25 




time [s] 

Figure 6: The temporal evolution of the temperature standard deviation for the ALT stirring protocol and 
flow configuration C of Table 2 for three modulation periods: 15 s, 30 s and 60 s. 

5.4- The evolution of extremum quantities 

The temporal evolutions of the maximum, mean and minimum temperature scalar values 
(and also some deduced values) are shown in Fig. 7 for the stirring flow C and the three 
stirring protocols: NM , CM and ALT. 

As mentioned in the section 2.4, we can verify (see Fig. 7(a)) that for each stirring pro- 
tocol, the mean fluid temperature evolves linearly with time with the same slope regardless 
of the flow conditions. Thus, the interesting quantities to follow with time are the maximum 
and minimum temperature values. We observe that after a transition period, the evolutions 
of the maximum and minimum temperatures also become linear, and evolve with the same 
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slope as the mean temperature regardless of the stirring protocol. The stirring protocol 
influences the gaps between the T^ aa , and T^ in curves and their respective position relative 
to the one. 

Hence, it is interesting to examine the temporal evolutions of T^ ax — T^ in plotted in 
Fig. 7(b), which shows that they have some shape similarity with those of a (see Fig. 4 (a)). 
In accordance with the results given by the temporal evolutions of the standard deviation of 
the temperature scalar (Fig. 4(a)), the ALT stirring protocol gives the minimum temperature 
difference T^ ax — T^ in and the NM stirring protocol gives the maximum one. Also, for the 
ALT stirring protocol T^ in has the closest values to T^. We recall that in the case of the 
fixed temperature boundary condition, the maximum temperature in the mixer is constant 
and only (Fig. 4(b) ) and the T^ in are increasing. 

In the case of heating and mixing with a fixed heat flux, the control of T^ in to prevent the 
persistence of cold spots is also problematic, and a particular effort has to be made to keep 
this minimal value as close as possible to the mean value. For this purpose, the last plot of 
Fig. 7 shows the temporal evolutions of the quantity T , ma ^ T » m . With this representation of 
the results, we can see the relative position of in the interval [T^ ax , T^ in ], and as before, 
the alternated stirring protocol always gives the best mixing and heat transfer efficiency For 
this case T^ ax — represents 80% of the total variation T^ ax — T^ in (T^ in is close to T^, 
see Fig. 7(a)). On the contrary, for the two other stirring protocols, T^ ax — is around 
50% of the total variation T^ ax - T* min . 

5.5. Probability distribution functions of the temperature scalars 

In this section, we examine the statistics of the temperature via the analysis of the 
probability distribution functions (PDFs) of the dimensionless temperature T* for the whole 
mixer section filled by the fluid. 

The PDFs of T* fields are shown in Fig. 8 for the flow configuration C, at different 
periodic times (for r = 30 s) and for the CM and ALT stirring protocols. In all of the 
PDFs, we notice the presence of a significant peak that corresponds to the most probable 
temperatures in the fluid. In the case of the imposed heat flux (Fig. 8(a)), while the PDFs of 
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Figure 7: For the flow configuration C of Table 2 and the three stirring protocols (NM, CM and ALT): (a) 
the temporal evolutions of T^ ax , and T^ in , (b) the temporal evolutions of T* lax — T^ in and (c) the 
temporal evolutions of (T^ ax - T^)/{T^ ax - T^ m ). 
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the CM stirring protocol are almost symmetrical, those of the ALT stirring protocol do not 
have a left tail. This indicates that the coldest temperatures have disappeared. Furthermore, 
the ALT protocol shows a higher peak (more exactly a bi-peak and at that particular time of 
the period) that translates towards the hot temperatures as time evolves. In both cases (CM 
and ALT), the PDFs recover their forms from period to period. Although the most probable 
temperature of the CM protocol is higher than that of the ALT protocol, its extended left 
tail manifests a lack of homogenization. 

When a constant temperature is imposed as the boundary condition (Fig. 8(b)), both 
the CM and ALT protocols do not have left tails. The peaks of the ALT protocol are located 
to the right of the CM peaks, towards the higher temperatures, but with lower amplitudes 
for the early times. For this type of boundary condition, when the peaks move towards the 
limit value T* = 1 their heights increase to reach the probability value of almost 1. On the 
contrary, for the imposed heat flux, the PDFs keep the same distribution while translating 
along the T* axis, a fact already shown for the evolutions of T^ ax — T^ in and T^ ax — T^. 

In Fig. 9, are shown for the ALT stirring protocol the PDFs of the rescaled dimensionless 
temperature: 

X = (12) 

a 

These PDFs are superimposed when plotted for different times during the mixing process, 
but at the same phase of the period. This is the signature of a strange eigenmode [12, 10], 
which is characterized by the production of persistent (self-similar) patterns in the flow. 
These patterns arise from a combination of stretching, folding and thermal diffusion. In 
fact, the temperature distribution evolves within the fluid during a period r, but it regains 
its previous form of one period earlier. 

5.6. Nusselt number evolutions 

The local heat transfer rate is quantified by the calculation of the mean Nusselt number, 
as indicated in section 4. The temporal evolutions of the Nusselt numbers are given sep- 
arately for the rods and tank in Fig. 10(a) and for all of the walls in Fig. 10(b). At the 
beginning of the process, the flow is not thermally or hydrodynamically developed, thus the 
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Nusselt number is higher. After a short time, the Nusselt number decreases and reaches an 
asymptotic state. This asymptotic state corresponds to a constant value for the CM stirring 
protocol and to a periodic state for the ALT stirring protocol (with a period identical to 
that of the wall velocity modulation). 

Due to this periodic wall displacement, the ALT stirring protocols present effectively 
large oscillations of the Nusselt number evolutions, and these oscillations are mainly given 
by the flow behavior around the rods [see Fig. 10(a)]. When summing the contributions of 
all the walls [Fig. 10(b)], the mean value of Nu given by the ALT stirring protocol is always 
greater than that given by the CM protocol (by a ratio of about 1.5), even if the walls have 
periodic stops in the ALT protocol. It is worth to note that the velocities of all the rotating 
walls are always greater in the case of CM protocol than of ALT one (Fig. 2), however, the 
Nusselt number is more important for the latter protocol: the parietal heat transfer is more 
governed by the patterns and the direction of the flow than by the magnitude of the velocity. 
Moreover, all of the flow configurations (A, B and C) give comparable results for the ALT 
protocol, while for the CM protocol, the C configuration is gives a lower Nu. 

5.7. Temperature gradient evolutions 

In Fig. 11, we present the evolution of the temperature gradients for the CM and ALT 
stirring protocols. The temperature scalar dissipation indicator \ g defined in section 4 is used 
for this purpose. In the case of the imposed heat flux (Fig. 11(a)), the temporal evolutions of 
X g are similar to those of the temperature standard deviation for the same stirring protocols 
(see Fig. 4(a)): after a first phase of gradient creation, its value stabilizes. However, the CM 
and ALT stirring protocols present large oscillations around their asymptotic value. The 
gradients of the temperature scalar do not vanish in time as is the case for mixing with a 
constant wall temperature boundary condition (Fig. 11(b)). In the latter case, the decrease of 
the temperature gradients is exponential, and again its rate is steeper for the ALT protocol. 
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6. Conclusion 

We have investigated, by numerical simulations, the coupled mixing and heating perfor- 
mances induced by chaotic advection in a 2D two-rod mixer for a constant wall heat flux 
boundary condition. Three different stirring protocols were compared: a non-modulated, a 
continuously modulated and an alternated (non-continuous) one for different flow configu- 
rations (depending on the respective direction of the rotation of the rod and tank walls). 
In order to study the thermal mixing enhancement mechanism within the fluids, different 
mixing indicators and statistical tools were used. According to the wall boundary condition 
considered (constant wall heat flux), the following main conclusions can be made based on 
the results obtained: 

• The most effective stirring protocol for thermal mixing is the alternated one (a con- 
clusion also obtained for the constant wall temperature boundary condition [12]). 

• For this alternated stirring protocol, the difference in thermal mixing induced by the 
choice of the flow configuration (A, B or C) is not very sensitive, contrary to the case 
of the continuously modulated stirring protocol. 

• The alternated stirring protocol seems to quickly suppress the coldest temperature 
during the mixing process, which is again not the case for the continuously modulated 
stirring protocol. 

• A great enhancement of heat transfer is thus obtained for the alternated stirring pro- 
tocol, confirmed by the Nusselt number evolutions. Despite the oscillations observed 
for the Nusselt number and temperature gradients evolutions, the homogeneity in the 
temperature during the thermal chaotic mixing is conserved (i.e. constant temperature 
standard deviation). 

We also showed the major differences between the two modes of parietal heating (fixed 
temperature or fixed heat flux) in the context of chaotic mixing. These differences can be 
summer ized in Tab. 3. 
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Fixed temperature 



Fixed heat flux 



Heat The amount of absorbed heat by the The amount of absorbed heat if fixed, 
absorbed fluid (and T m ) depends on the flow. 

T m The evolution of T m is asymptotic. The evolution of T m is linear and con- 



tinuously evolving. 



A high degree of homogenization 
(very low a) can be achieved as the 
mixing process is continued, a de- 
creases exponentially with time. 



The degree of homogenization settles 
at a fixed value. 



Asymtotic For T m or T* 
thermal 
regime 



For a, T* max 



T* ■ or T* 

mm max 



-^m> Xg- 



Mixing The efficiency of the mixing protocol 
efficiency impacts the rate of homogenization 
and the decay rate of the tempera- 
ture gradient \g- 



The efficiency of the mixing protocol 
impacts the difference between the 
extrema of the temperatures, as well 
as a and \g- All these values remain 
constant as the mixing proceeds. 



Table 3: Major differences between the fixed temperature and the fixed heat flux modes of parietal heating 
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The mixing and heating using a prescribed heat flux always gives a non-homogenized 
fluid temperature. On the other hand, heating with a fixed temperature suffers from the 
asymptotic evolution of the quantity of heat transmitted to the fluid. In an industrial process, 
the heating of a highly viscous fluid can be first started with a fixed flux to benefit from a 
high rate of energy supply. An alternated protocol will ensure a low gradient in the fluid and 
will help to maintain the maximal temperature below the fluid damage temperature. When 
the target temperature is approached, it can be fixed at the mixer walls to finish the mixing 
with the suitable degree of homogenization. 

Further studies will explore the effect of the temperature-dependence of the physical 
properties on the efficiency of the thermal chaotic mixing. 
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Appendix A: The numerical method 

To describe the discretization practice used in the code, we can write the above equations 
(3 and 6) in the generic convection-diffusion form: 



d_ 

Of 



[ p<pdV+ [ P <pU-ndS= [ rV0-ndS+ I S^dV (13) 
Jv Js Js Jv 



where is a source term. The spatial schemes approximating the diffusive and convective 
fluxes are both second-order accurate. The discretization of the diffusion term is performed 
by approximating the surface integrals by taking the sum over all the cell faces / (see Figure 
12): 

[rV<f ) -ndS = Y, TfA f (vfof ■ ftf (14) 
Js f 

In unstructured meshes, non-orthogonality is not an exception and it needs to be correctly 
handled. Thus, the normal gradient (V0)/- • n/ is decomposed into an implicit contribution 
that uses the values of <p at the centers of the two cells sharing the face / (the first term 
in the right hand side (RHS) of Eq. (15) ) and a non-orthogonality correction term treated 
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explicitly by a deferred approach in order to preserve the second-order accuracy of this 
centered differencing. We use here the over-relaxed decomposition as suggested by Jasak 
[24] to enhance the convergence properties of the diffusive term discretization: 



(V*)/ • % = ^ a r c ?\ + V* • I n f - -=^- ] (15) 
\\d\\ d-n } \ d-n f/ 

dis the vector joining the centers of the two cells (see Figure 12). The average gradient V0 
is interpolated from the gradients of these neighboring cells. 

The gradients of the variables at the cell centers are computed by Gauss' theorem: 

/ 

where 0/ is the mean value of the variable using the values at the centers of two cells sharing 
the face /: 

/ = £0 c + (l-O<K with £== (17) 

Once the gradient is calculated at all the computational cells, its new values are used to give 
a new estimate of <bf as: 



d 



(<K + v$ Nb ■ iv 6 /) + + v$ c ■ cf) 

These new values of 0/ are used to recompute the gradients more accurately [25] . 

Convection terms are also transformed to a sum over the faces / composing the surface 

S: 

[ p( f>U-ndS = J2(pM)fU r n f (19) 
Js f 

where the face values <pf require an appropriate interpolation to achieve accuracy and bound- 
edness. For this purpose, we use the non-linear high resolution (HR) bounded scheme CU- 
B1STA of Alves et al. [26], in the 7 formulation of Ng et al. [27], where they express <frf as a 
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function of the upwind value of <f) and of its centered differencing (CD) value: 



HR 

7 




(20) 



The coefficient 7 is determined face by face based on the local shape of the flow solution, using 
the normalized variable diagram (NVD) framework and observing the convection bounded- 
ness criterion (CBC)[28]. The first term of the RHS of Eq. (20) is accounted for implicitly, 
while the second term is treated explicitly with the deferred-correction practice. 

The pressure- velocity coupling is ensured by the SIMPLE algorithm [29], while the mass 
fluxes at the cell faces are evaluated using the Rhie-Chow interpolation [30] to avoid pressure 
checkerboarding. The implicit three-time step Gear's scheme of second-order accuracy is used 
to discretize the unsteady terms: 



The superscript n stands for the current time step and At = t n — t n 1 for the time step. The 
RHS of Eq. (13) is taken at time t n . 

At each iteration, the discretization technique presented above leads to a linear system 
with a non-symmetric sparse matrix for each variable. These linear systems are solved using 
an ILU-preconditioned GMRES solver using the implementation of the IML++ library [31]. 

Appendix B: The code validation 

Different validation test cases have been performed to assess the accuracy of the code. 
We present here two cases corresponding to enclosed flows with sliding boundaries that make 
them similar to the flow in our two-rod mixer. They are both 2D flows. The first case is 
the classical lid-driven square cavity flow at a Reynolds number of Re = 1000. Our findings 
are compared to the numerical results of Botella and Peyret [32] , which were obtained using 
a spectral Chebyshev method, and to those of Ghia et al. [33], who used a vorticity-stream 
function formulation. The second test case has a circular sliding boundary it is the polar 
lid-driven cavity flow studied experimentally and numerically by Fuchs and Tillmark [34]. 




(21) 
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We compare here our results to their experimental measurements obtained by laser Doppler 
anemometry at Reynolds numbers Re = 60 and 350. 

The meshes used in the two cases are shown in Fig. 13. They are two-dimensional (i.e. 
composed of one layer of cells) and of two different types. The mesh of the square cavity is 
composed of triangular-based prisms, only to demonstrate the ability of the code to handle 
non-orthogonal grids. It contains 10468 computational cells and is refined near the walls. 
The mesh of the polar cavity benefits from the hybrid mesh support and uses quadrilateral- 
based prisms (hexahedrals) near the boundaries to increase the mesh resolution and enhance 
the orthogonality, while triangular-based prisms are used elsewhere to ensure a flexible and 
homogeneous distribution of cells. This mesh is composed of 9610 cells. 

In Fig. 14 we present the obtained results and their comparisons. For the square cavity 
(Fig. 14(a)), the x velocity component u is plotted along the vertical central axis and the y 
velocity component v is plotted along the horizontal central axis. Our results match quite 
well with the benchmark results of Botella and Peyret [32] and Ghia et al. [33]. Fig. 14(b) 
shows the angular velocity along the central radial axis of the polar cavity. Good agreement 
is obtained with the experimental results of Fuchs and Tillmark [34]. Supplementary details 
and validation tests of this code are related in [21, 22]. 
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Figure 8: PDFs of the dimensionlcss temperature T* at different periodic times for the CM and ALT stirring 
protocols and flow configuration C (r = 30 s). (a) The fixed heat flux boundary condition, and (b) the fixed 
temperature boundary condition. 
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Figure 9: PDFs of the rescaled dimensionlcss temperature X at different periodic times for the ALT stirring 
protocols and the flow configuration C (r = 30 s). 
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Figure 10: The Nusselt number evolutions: (a) for the NM, CM and ALT stirring protocols with the rod 
and tank contributions separated (flow configuration C), and (b) for the flow configurations A, B and C and 
for the CM and ALT stirring protocols with the rod and tank contributions combined. 
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Figure 11: The temporal evolutions of the scalar temperature dissipation indicator for the NM, CM and 
ALT stirring protocols with flow configuration C. (a) The imposed heat flux, (b) The imposed temperature 
(logarithmic scale). 
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Figure 12: A computational cell C and one of its neighbors JV& 




Figure 13: Computational meshes used for the validation, (a) A lid-driven square cavity and (b) a polar 
cavity with a rotating inner circular wall. 
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Figure 14: (a) The lid-driven square cavity: the u and v velocity components along the vertical and horizontal 
axes, Re = 1000. (b) The polar cavity: the ug velocity component along radial axis, Re = 60 and 350. 
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